Groundwater discharge
This tutorial is very similar to the river discharge notebook
In this tutorial we will simulate a fictitious radioactive tracer that is injected into the ocean by groundwater discharge. The 200 major rivers dataset from Luijendijk et al. (2020) is available from within the AIBECS. Once "born", our ficitious tracer decays with a parameter timescale $\tau$ as it flows through ocean basins.
The 3D tracer equation is:
where $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \nabla \right]$ represents the ocean circulation transport. (Tracer transport operators are described in the documentation.) The riverine source of the tracer is $s_\mathsf{gw}$, and $x / \tau$ is the decay rate.
In AIBECS, we must recast this equation in the generic form
We start by telling Julia that we want to use the AIBECS and the OCIM2 transport matrix for the ocean circulation.
using AIBECS
grd, T_OCIM2 = OCIM2.load()(,
[1 , 1] = 0.000197784
[2 , 1] = 2.34279e-9
[10384 , 1] = -1.95995e-7
[10442 , 1] = -0.000191612
[10443 , 1] = 4.80961e-9
[20825 , 1] = -1.83059e-9
[20883 , 1] = 5.00768e-9
[1 , 2] = -5.02516e-8
[2 , 2] = 0.000187531
⋮
[200160, 200159] = -2.19656e-8
[197886, 200160] = 1.08199e-10
[199766, 200160] = 6.70981e-9
[199777, 200160] = -1.26352e-9
[199778, 200160] = -3.39279e-9
[199779, 200160] = 7.59316e-9
[199790, 200160] = -7.41018e-9
[200156, 200160] = -3.44106e-8
[200159, 200160] = -2.00303e-8
[200160, 200160] = 5.27945e-8, Quantity{Float64,𝐍 𝐋⁻³ 𝐓⁻¹,Unitful.FreeUnits{(m⁻³, mol, yr⁻¹),𝐍 𝐋⁻³ 𝐓⁻¹,nothing}}[1.82368935574118e-11 mol m⁻³ yr⁻¹, 1.9469772194929453e-11 mol m⁻³ yr⁻¹, 3.4879359939844136e-13 mol m⁻³ yr⁻¹, 1.5360041759458074e-12 mol m⁻³ yr⁻¹, 1.47673234957578e-12 mol m⁻³ yr⁻¹, 5.300237184679614e-11 mol m⁻³ yr⁻¹, 8.595855761423191e-11 mol m⁻³ yr⁻¹, 1.4013789381514675e-10 mol m⁻³ yr⁻¹, 2.0592425177423245e-10 mol m⁻³ yr⁻¹, 2.3948293487252974e-10 mol m⁻³ yr⁻¹ … 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹], Quantity{Float64,𝐍 𝐋⁻³ 𝐓⁻¹,Unitful.FreeUnits{(m⁻³, mol, yr⁻¹),𝐍 𝐋⁻³ 𝐓⁻¹,nothing}}[1.2487363516886739e-5 mol m⁻³ yr⁻¹, 1.333155354686199e-5 mol m⁻³ yr⁻¹, 2.388297356860729e-7 mol m⁻³ yr⁻¹, 1.0517494357308457e-6 mol m⁻³ yr⁻¹, 1.0111641880370838e-6 mol m⁻³ yr⁻¹, 3.629235880686242e-5 mol m⁻³ yr⁻¹, 5.8858475701302695e-5 mol m⁻³ yr⁻¹, 9.595673830367888e-5 mol m⁻³ yr⁻¹, 0.0001410026867104604 mol m⁻³ yr⁻¹, 0.00016398135211074882 mol m⁻³ yr⁻¹ … 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹, 0.0 mol m⁻³ yr⁻¹])The transport is
T_gw(p) = T_OCIM2T_gw (generic function with 1 method)
For the radioactive decay, we simply use
function decay(x, p)
@unpack τ = p
return x / τ
enddecay (generic function with 1 method)
To build the river sources, we will load the geographic locations and discharge of groundwater (in m³ yr⁻¹) from the Luijendijk et al. (2020)) dataset.
gws = GroundWaters.load()40409-element Array{AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}},1}:
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(35.1794887415, 24.2257701314, 47358.8148595 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(35.3417438136, 23.8746131396, 185142.093516 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(35.4161096656, 23.6050233764, 171000.541515 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(35.4266133609, 23.4030494462, 206889.096486 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(35.7820677063, 22.6137148344, 667545.345579 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(35.8962644211, 22.6987844211, 72625.201875 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(35.8614759879, 22.5822032466, 911480.517723 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(36.6601448123, 21.8794834664, 472586.358228 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(37.2143679963, 19.334562086, 1.33262583591e6 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(37.3177380952, 19.3185952381, 28384.2888 m³ yr⁻¹)
⋮
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(144.157841392, 72.2055303019, 11.7123541475 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(144.236000631, 72.2728170091, 41.4176247705 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(145.145353536, 71.8562156581, 73.1538069047 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(128.065238136, 73.2051409179, 105.676353367 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(143.155736084, 75.5195182604, 93.4943801401 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(143.243027112, 75.6684283401, 48.7752847061 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(142.916389905, 75.6972027214, 69.4594679685 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(128.501004892, 73.1684350529, 87.5895493016 m³ yr⁻¹)
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(141.832103547, 75.6537345857, 45.9949603817 m³ yr⁻¹)This is an array of groundwater sources, for which the type GroundWaterSource{T} contains the lat–lon coordinates and discharge in m³ yr⁻¹. For example, the first element is
gws[1]AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64,𝐋³ 𝐓⁻¹,Unitful.FreeUnits{(m³, yr⁻¹),𝐋³ 𝐓⁻¹,nothing}}}(35.1794887415, 24.2257701314, 47358.8148595 m³ yr⁻¹)We can check the locations with
using Plots
scatter([x.lon for x in gws], [x.lat for x in gws],
zcolor=log10.(ustrip.([x.VFR for x in gws] / u"m^3/s")),
clim=(0,7), colorbartitle="log₁₀(discharge / (1 m³ s⁻¹))")
We can regrid these into the OCIM2 grid and return the corresponding vector with
gws_grd = ustrip.(u"m^3/s", regrid(gws, grd))200160-element Array{Float64,1}:
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
⋮
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0(We convert years to seconds for AIBECS)
We can control the global magnitude of the discharge by chosing a concentration of our tracer in groundwater, $C_\mathsf{gw}$, in mol m⁻³. For that, we need the volumes of the grid boxes
v = volumevec(grd)200160-element Array{Float64,1}:
5.693492696598497e11
6.267656666812283e11
6.834351352130837e11
7.392901411788552e11
7.942641211123751e11
8.482915614828765e11
9.013080767686609e11
9.532504861865775e11
1.004056888985702e12
1.053666738215696e12
⋮
2.7437236229729617e13
2.6923849140323855e13
2.83653451428806e13
2.7917925846214992e13
2.7437236229729617e13
2.6923849140323855e13
2.83653451428806e13
2.7917925846214992e13
2.7437236229729617e13because the source is given by the mass (discharge × concentration) divided by the volume of each box:
function s_gw(p)
@unpack C_gw = p
return @. gws_grd * C_gw / v
ends_gw (generic function with 1 method)
We then write the generic $\boldsymbol{G}$ function, which is
G_gw(x,p) = s_gw(p) - decay(x,p)G_gw (generic function with 1 method)
Parameters
We specify some initial values for the parameters and also include units.
import AIBECS: @units, units
import AIBECS: @initial_value, initial_value
@initial_value @units struct GroundWatersParameters{U} <: AbstractParameters{U}
τ::U | 20.0 | u"yr"
C_gw::U | 1.0 | u"mol/m^3"
endinitial_value (generic function with 53 methods)
Finally, thanks to the initial values we provided, we can instantiate the parameter vector succinctly as
p = GroundWatersParameters()│ Row │ Symbol │ Value │ Initial value │ Unit │ │ │ Symbol │ Float64 │ Float64 │ FreeUni… │ ├─────┼────────┼─────────┼───────────────┼──────────┤ │ 1 │ τ │ 20.0 │ 20.0 │ yr │ │ 2 │ C_gw │ 1.0 │ 1.0 │ mol m⁻³ │
We generate the state function F and its Jacobian ∇ₓF,
F, ∇ₓF = state_function_and_Jacobian(T_gw, G_gw)(AIBECS.var"#F#35"{typeof(Main.ex-6_groundwater_discharge.T_gw),typeof(Main.ex-6_groundwater_discharge.G_gw)}(Main.ex-6_groundwater_discharge.T_gw, Main.ex-6_groundwater_discharge.G_gw), AIBECS.var"#∇ₓF#37"{typeof(Main.ex-6_groundwater_discharge.T_gw),AIBECS.var"#∇ₓG#36"{typeof(Main.ex-6_groundwater_discharge.G_gw)}}(Main.ex-6_groundwater_discharge.T_gw, AIBECS.var"#∇ₓG#36"{typeof(Main.ex-6_groundwater_discharge.G_gw)}(Main.ex-6_groundwater_discharge.G_gw)))generate the steady-state problem prob,
nb = sum(iswet(grd))
x = ones(nb) # initial guess
prob = SteadyStateProblem(F, ∇ₓF, x, p)SteadyStateProblem with uType Array{Float64,1}
u0: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]and solve it
s = solve(prob, CTKAlg()).u * u"mol/m^3"200160-element Array{Quantity{Float64,𝐍 𝐋⁻³,Unitful.FreeUnits{(m⁻³, mol),𝐍 𝐋⁻³,nothing}},1}:
-1.1976569114251033e-7 mol m⁻³
-1.447973438008902e-7 mol m⁻³
-1.4879453161968372e-7 mol m⁻³
-1.5817358327494972e-7 mol m⁻³
-1.952022620217737e-7 mol m⁻³
-2.0677685877025737e-7 mol m⁻³
-2.244366103876224e-7 mol m⁻³
-4.1123842884193485e-7 mol m⁻³
-5.599692326820449e-7 mol m⁻³
-4.425251364361234e-7 mol m⁻³
⋮
5.356179456771712e-8 mol m⁻³
9.890403277283724e-8 mol m⁻³
-1.061571573408199e-7 mol m⁻³
-3.777197030624626e-8 mol m⁻³
7.586448212525718e-8 mol m⁻³
1.8682291679180843e-7 mol m⁻³
-6.145591682153108e-8 mol m⁻³
-6.283597352127268e-8 mol m⁻³
6.068167250823727e-8 mol m⁻³Let's now run some visualizations using the plot recipes. Taking a horizontal slice of the 3D field at 200m gives
cmap = :viridis
plothorizontalslice(s, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,100))
and at 500m:
plothorizontalslice(s, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,25))
Or we can increase the decay timescale (×10) and decrease the groundwater concentration (÷10) to get a different (more well-mixed) tracer distribution:
p = GroundWatersParameters(τ = 200.0u"yr", C_gw = 0.1u"mol/m^3")
prob = SteadyStateProblem(F, ∇ₓF, x, p)
s_τ50 = solve(prob, CTKAlg()).u * u"mol/m^3"
plothorizontalslice(s_τ50, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,25))
This page was generated using Literate.jl.